/************************************************************************************************
*			Primeiro Trabalho de Programacao Concorrente				*	
*	Implementacao paralela do algoritmo de Gauss-Legendre para calculo de casas do pi	*
*												*
*	Alunos:											*
*		Felipe Augusto do Amaral (7239221)						*
*		Julio Cesar Gaban (7152872)							*
*		Tiago Santana de Nazare (7153118)						*
*												*
*	Execucao:										*
*		./<executavel> <numero_de_iteracoes>						*
*												*
*	Compilacao:										*
*		gcc <fonte> -o <executavel> -lgmp -pthread -ansi -Wall				*
*************************************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <pthread.h>

#define PRECISAO 33219281

/* Parametros para o calculo de ak */
struct calc_a_t{
	mpf_t *a1;
	mpf_t *b1;
	mpf_t *a2;
};
typedef struct calc_a_t calc_a_t;

/* Parametros para o calculo de bk */
struct calc_b_t{
	mpf_t *a1;
	mpf_t *b1;
	mpf_t *b2;
};
typedef struct calc_b_t calc_b_t;

/* Parametros para o calculo de tk */
struct calc_t_t{
	mpf_t *t1;
	mpf_t *p1;
	mpf_t *a1;
	mpf_t *a2;
	mpf_t *t2;
};
typedef struct calc_t_t calc_t_t;

/* Parametros para o calculo de pk */
struct calc_p_t{
	mpf_t *p1;
	mpf_t *p2;
};
typedef struct calc_p_t calc_p_t;

void *calc_a(void *param);
void *calc_b(void *param);
void *calc_t(void *param);
void *calc_p(void *param);

int main(int argc, char *argv[]){
	unsigned long int i, nIteracoes;
	mpf_t a1, a2, b1, b2, p1, p2, t1, t2, pi;
	pthread_t thread_calc_a, thread_calc_b, thread_calc_t, thread_calc_p;
	calc_a_t data_calc_a;
	calc_b_t data_calc_b;
	calc_p_t data_calc_p;
	calc_t_t data_calc_t;

	/* Verifica se o argumento de entrada (numero de iteracoes) foi passado */
	if(argc != 2){
		printf("Argumentos invalidos!!!\n");
		printf("Use: %s <numero_de_iteracoes>\n", argv[0]);
		exit(1);
	}

	nIteracoes = strtoul(argv[1], NULL, 0);

	mpf_set_default_prec(PRECISAO);

	/* Inicia as variaveis */
	/* a0 = 1 */
	mpf_init_set_d(a1, 1.0);

	/* b0 = 1 / 2^(0.5) */
	mpf_init_set_d(b1, 2.0);
	mpf_sqrt(b1, b1);
	mpf_div(b1, a1, b1);

	/* t0 = 1 / 4 */
	mpf_init_set_d(t1, 0.25);

	/* p0 = 1 */
	mpf_init_set_d(p1, 1.0);

	mpf_init(a2);
	mpf_init(b2);
	mpf_init(t2);
	mpf_init(p2);
	mpf_init(pi);

	data_calc_a.a1 = &a1;
	data_calc_a.b1 = &b1;
	data_calc_a.a2 = &a2;

	data_calc_b.a1 = &a1;
	data_calc_b.b1 = &b1;
	data_calc_b.b2 = &b2;

	data_calc_t.t1 = &t1;
	data_calc_t.p1 = &p1;
	data_calc_t.a1 = &a1;
	data_calc_t.a2 = &a2;
	data_calc_t.t2 = &t2;

	data_calc_p.p1 = &p1;
	data_calc_p.p2 = &p2;

	/* Faz as iteracoes do algoritmo */
	for(i = 0; i < nIteracoes; ++i){
		/* Cria as threads para calcular ak+1, bk+1, pk+1 */
		pthread_create(&thread_calc_a, NULL, calc_a, (void *) &data_calc_a);
		pthread_create(&thread_calc_b, NULL, calc_b, (void *) &data_calc_b);
		pthread_create(&thread_calc_p, NULL, calc_p, (void *) &data_calc_p);

		/* Espera ak ser calculado para calcular tk+1 (tk+1 precisa de ak+1) */
		pthread_join(thread_calc_a, NULL);
		pthread_create(&thread_calc_t, NULL, calc_t, (void *) &data_calc_t);

		pthread_join(thread_calc_b, NULL);
		pthread_join(thread_calc_p, NULL);
		pthread_join(thread_calc_t, NULL);

		mpf_set(a1, a2);
		mpf_set(b1, b2);
		mpf_set(p1, p2);
		mpf_set(t1, t2);
	}

	/* pi = ((a2 + b2)^2) / 4t2 */
	mpf_add(pi, a2, b2);
	mpf_pow_ui(pi, pi, 2);
	mpf_div_ui(pi, pi, 4);
	mpf_div(pi, pi, t2);

	mpf_out_str(stdout, 10, 0, pi);

	mpf_clear(a1);
	mpf_clear(a2);	
	mpf_clear(b1);
	mpf_clear(b2);
	mpf_clear(p1);
	mpf_clear(p2);
	mpf_clear(pi);

	return 0;
}

/* Calcula ak+1 */
void *calc_a(void *param){
	calc_a_t data = *((calc_a_t *)param);

	/* a2 = (a1 + b1) / 2 */
	mpf_add(*data.a2, *data.a1, *data.b1);
	mpf_div_ui(*data.a2, *data.a2, 2);

	pthread_exit(NULL);
}

/* Calcula bk+1 */
void *calc_b(void *param){
	calc_b_t data = *((calc_b_t *)param);

	/* b2 = (a1 * b1)^(0.5) */
	mpf_mul(*data.b2, *data.a1, *data.b1);
	mpf_sqrt(*data.b2, *data.b2);

	pthread_exit(NULL);
}

/* Calcula tk+1 */
void *calc_t(void *param){
	calc_t_t data = *((calc_t_t *)param);

	/* t2 = t1 - p1((a1 - a2)^2) */
	mpf_sub(*data.t2, *data.a1, *data.a2);
	mpf_pow_ui(*data.t2, *data.t2, 2);
	mpf_mul(*data.t2, *data.p1, *data.t2);
	mpf_sub(*data.t2, *data.t1, *data.t2);

	pthread_exit(NULL);
}

/* Calcula pk+1 */
void *calc_p(void *param){
	calc_p_t data = *((calc_p_t *)param);

	/* p2 = 2 * p1 */
	mpf_mul_ui(*data.p2, *data.p1, 2);

	pthread_exit(NULL);
}

